############
# packages #
############
library(PanelMatch)
library(plm)
library(dplyr)
library(tidyr)
library(data.table)
library(ggplot2)
library(patchwork)

###############
# set-up data #
###############
# read data set
dat.anlsys <- read.csv("Party data for analysis.csv")
dat.anlsys <- as.data.frame(dat.anlsys)
# turn variables from integer to numeric
dat.anlsys <- dat.anlsys %>% mutate_if(is.integer, as.numeric)
# keep time and unit id as integer
dat.anlsys$period <- as.integer(dat.anlsys$period)
dat.anlsys$v2paid <- as.integer(dat.anlsys$v2paid)

###################################################################
# select matching refinement method - moderate reform (Figure A1) #
###################################################################
## generate panel data for PanelMatch
dat.anlsys.dlgts <- PanelData(panel.data = dat.anlsys, time.id = "period",
                            unit.id = "v2paid",outcome = "v2pavote",
                            treatment = "prmrs.dlgts")
## match sets based on different matching methods
# no method (equal weights)
ref.none.dlgts <- PanelMatch(panel.data = dat.anlsys.dlgts, lag = 2, match.missing = TRUE, 
                             refinement.method = "none",
                             covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                               v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                               I(lag(gov.mmbr,1)) + I(lag(prmrs.dlgts.swtch.other,1)) + 
                               e_gdppc + v2paallian, 
                             qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# mahalanobis distance
ref.mal.dlgts <- PanelMatch(panel.data = dat.anlsys.dlgts, lag = 2, match.missing = TRUE, 
                            refinement.method = "mahalanobis",
                            covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                              v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                              I(lag(gov.mmbr,1)) + I(lag(prmrs.dlgts.swtch.other,1)) + 
                              e_gdppc + v2paallian, 
                            qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# Propensity score matching
ref.psm.dlgts <- PanelMatch(panel.data = dat.anlsys.dlgts, lag = 2, match.missing = TRUE, 
                            refinement.method = "ps.match",
                            covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                              v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                              I(lag(gov.mmbr,1)) + I(lag(prmrs.dlgts.swtch.other,1)) + 
                              e_gdppc + v2paallian, 
                            qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# Propensity score weighing
ref.psw.dlgts <- PanelMatch(panel.data = dat.anlsys.dlgts, lag = 2, match.missing = TRUE, 
                            refinement.method = "ps.weight",
                            covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                              v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                              I(lag(gov.mmbr,1)) + I(lag(prmrs.dlgts.swtch.other,1)) + 
                              e_gdppc + v2paallian, 
                            qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)

## compare balance stats
ref.none.dlgts.blnce <- get_covariate_balance(ref.none.dlgts, panel.data = dat.anlsys.dlgts,
                      covariates = c("v2pavote", "v2elparlel","v2pariglef_ord","v2pawomlab_ord","v2paimmig_ord",
                                     "e_gdppc", "v2paallian","prmrs.dlgts.swtch.other","gov.mmbr"))
ref.mal.dlgts.blnce <- get_covariate_balance(ref.mal.dlgts, panel.data = dat.anlsys.dlgts,
                      covariates = c("v2pavote", "v2elparlel","v2pariglef_ord","v2pawomlab_ord","v2paimmig_ord",
                                     "e_gdppc", "v2paallian","prmrs.dlgts.swtch.other","gov.mmbr"))
ref.psm.dlgts.blnce <- get_covariate_balance(ref.psm.dlgts, panel.data = dat.anlsys.dlgts,
                      covariates = c("v2pavote", "v2elparlel","v2pariglef_ord","v2pawomlab_ord","v2paimmig_ord",
                                     "e_gdppc", "v2paallian","prmrs.dlgts.swtch.other","gov.mmbr"))
ref.psw.dlgts.blnce <- get_covariate_balance(ref.psw.dlgts, panel.data = dat.anlsys.dlgts,
                      covariates = c("v2pavote", "v2elparlel","v2pariglef_ord","v2pawomlab_ord","v2paimmig_ord",
                                     "e_gdppc", "v2paallian","prmrs.dlgts.swtch.other","gov.mmbr"))
print(c(ref.none.dlgts.blnce, ref.mal.dlgts.blnce, ref.psm.dlgts.blnce, ref.psw.dlgts.blnce))

## plot balance before and after refinement (by each method)
# generate dataframe for plot + reshape for plot
ref.none.dlgts.blnce.df <- as.data.frame(ref.none.dlgts.blnce)
ref.none.dlgts.blnce.df <- setDT(ref.none.dlgts.blnce.df, keep.rownames = TRUE)[]
ref.none.dlgts.blnce.long <- gather(ref.none.dlgts.blnce.df, "covariate", "balance",
                                    ref.none.dlgts.att.v2pavote:ref.none.dlgts.att.gov.mmbr)
ref.none.dlgts.blnce.long$rn <- factor(ref.none.dlgts.blnce.long$rn, levels=c("t_2","t_1","t_0"))
ref.none.dlgts.blnce.long$dv <- ifelse(
  ref.none.dlgts.blnce.long$covariate=="ref.none.dlgts.att.v2pavote",1,0)
ref.none.dlgts.blnce.long <- filter(ref.none.dlgts.blnce.long,rn!="t_0")

ref.mal.dlgts.blnce.df <- as.data.frame(ref.mal.dlgts.blnce)
ref.mal.dlgts.blnce.df <- setDT(ref.mal.dlgts.blnce.df, keep.rownames = TRUE)[]
ref.mal.dlgts.blnce.long <- gather(ref.mal.dlgts.blnce.df, "covariate", "balance",
                                   ref.mal.dlgts.att.v2pavote:ref.mal.dlgts.att.gov.mmbr)
ref.mal.dlgts.blnce.long$rn <- factor(ref.mal.dlgts.blnce.long$rn, levels=c("t_2","t_1","t_0"))
ref.mal.dlgts.blnce.long$dv <- ifelse(
  ref.mal.dlgts.blnce.long$covariate=="ref.mal.dlgts.att.v2pavote",1,0)
ref.mal.dlgts.blnce.long <- filter(ref.mal.dlgts.blnce.long,rn!="t_0")


ref.psm.dlgts.blnce.df <- as.data.frame(ref.psm.dlgts.blnce)
ref.psm.dlgts.blnce.df <- setDT(ref.psm.dlgts.blnce.df, keep.rownames = TRUE)[]
ref.psm.dlgts.blnce.long <- gather(ref.psm.dlgts.blnce.df, "covariate", "balance",
                                   ref.psm.dlgts.att.v2pavote:ref.psm.dlgts.att.gov.mmbr)
ref.psm.dlgts.blnce.long$rn <- factor(ref.psm.dlgts.blnce.long$rn, levels=c("t_2","t_1","t_0"))
ref.psm.dlgts.blnce.long$dv <- ifelse(
  ref.psm.dlgts.blnce.long$covariate=="ref.psm.dlgts.att.v2pavote",1,0)
ref.psm.dlgts.blnce.long <- filter(ref.psm.dlgts.blnce.long,rn!="t_0")

ref.psw.dlgts.blnce.df <- as.data.frame(ref.psw.dlgts.blnce)
ref.psw.dlgts.blnce.df <- setDT(ref.psw.dlgts.blnce.df, keep.rownames = TRUE)[]
ref.psw.dlgts.blnce.long <- gather(ref.psw.dlgts.blnce.df, "covariate", "balance",
                                   ref.psw.dlgts.att.v2pavote:ref.psw.dlgts.att.gov.mmbr)
ref.psw.dlgts.blnce.long$rn <- factor(ref.psw.dlgts.blnce.long$rn, levels=c("t_2","t_1","t_0"))
ref.psw.dlgts.blnce.long$dv <- ifelse(
  ref.psw.dlgts.blnce.long$covariate=="ref.psw.dlgts.att.v2pavote",1,0)
ref.psw.dlgts.blnce.long <- filter(ref.psw.dlgts.blnce.long,rn!="t_0")

# plot balance (Figure A1)
ref.none.dlgts.blnce.plot <- ggplot(dat=ref.none.dlgts.blnce.long, aes(x=rn, y=balance, alpha=as.factor(dv),
                                                                       group=covariate)) +
  geom_line(aes(size=as.factor(dv))) + 
  geom_hline(yintercept=0.25, linetype="dashed") +
  geom_hline(yintercept=-0.25, linetype="dashed") +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_x_discrete(labels=c("t-2", "t-1")) +
  scale_alpha_discrete(range = c(0.3, 1), labels=c("Covariates","Vote Share")) +
  scale_size_discrete(range = c(0.5,1), labels=c("Covariates","Vote Share")) + 
  ylim(-1.6,1.6) + 
  labs(title="No Refinement", x="Elections Before Treatment", y="Standardzies Mean \nDifferences") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5),
  )

ref.mal.dlgts.blnce.plot <- ggplot(dat=ref.mal.dlgts.blnce.long, aes(x=rn, y=balance, alpha=as.factor(dv),
                                                                     group=covariate)) +
  geom_line(aes(size=as.factor(dv))) + 
  geom_hline(yintercept=0.25, linetype="dashed") +
  geom_hline(yintercept=-0.25, linetype="dashed") +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_x_discrete(labels=c("t-2", "t-1")) +
  scale_alpha_discrete(range = c(0.3, 1), labels=c("Covariates","Vote Share")) +
  scale_size_discrete(range = c(0.5,1), labels=c("Covariates","Vote Share")) + 
  ylim(-1.6,1.6) + 
  labs(title="Mahalanobis", x="Elections Before Treatment", y="Standardzies Mean \nDifferences") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5),
  )

ref.psm.dlgts.blnce.plot <- ggplot(dat=ref.psm.dlgts.blnce.long, aes(x=rn, y=balance, alpha=as.factor(dv),
                                                                     group=covariate)) +
  geom_line(aes(size=as.factor(dv))) + 
  geom_hline(yintercept=0.25, linetype="dashed") +
  geom_hline(yintercept=-0.25, linetype="dashed") +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_x_discrete(labels=c("t-2", "t-1")) +
  scale_alpha_discrete(range = c(0.3, 1), labels=c("Covariates","Vote Share")) +
  scale_size_discrete(range = c(0.5,1), labels=c("Covariates","Vote Share")) + 
  ylim(-1.6,1.6) + 
  labs(title="Propensity Score\nMatching", x="Elections Before Treatment", y="Standardzies Mean \nDifferences") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5),
  )

ref.psw.dlgts.blnce.plot <- ggplot(dat=ref.psw.dlgts.blnce.long, aes(x=rn, y=balance, alpha=as.factor(dv),
                                                                     group=covariate)) +
  geom_line(aes(size=as.factor(dv))) + 
  geom_hline(yintercept=0.25, linetype="dashed") +
  geom_hline(yintercept=-0.25, linetype="dashed") +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_x_discrete(labels=c("t-2", "t-1")) +
  scale_alpha_discrete(range = c(0.3, 1), labels=c("Covariates","Vote Share")) +
  scale_size_discrete(range = c(0.5,1), labels=c("Covariates","Vote Share")) + 
  ylim(-1.6,1.6) + 
  labs(title="Propensity Score\nWeighting", x="Elections Before Treatment", y="Standardzies Mean \nDifferences") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5),
  )

plot.dlgts.blnce <- ref.none.dlgts.blnce.plot + ref.mal.dlgts.blnce.plot + ref.psm.dlgts.blnce.plot + ref.psw.dlgts.blnce.plot +  
  plot_layout(ncol = 2, guides = "collect") & theme(legend.position = 'bottom')
# save plot
ggsave(file="FigureA1.pdf", plot=plot.dlgts.blnce, width=7, height=7)

#################################################################
# select matching refinement method - strict reform (Figure A2) #
#################################################################
# generate panel data for PanelMatch
dat.anlsys.mmbrs <- PanelData(panel.data = dat.anlsys, time.id = "period",
                              unit.id = "v2paid",outcome = "v2pavote",
                              treatment = "prmrs.mmbrs")
## match sets based on different matching methods
# no method (equal weights)
ref.none.mmbrs <- PanelMatch(panel.data = dat.anlsys.mmbrs, lag = 2, match.missing = TRUE, 
                             refinement.method = "none",
                             covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                               v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                               I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                               e_gdppc + v2paallian, 
                             qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# mahalanobis distance
ref.mal.mmbrs <- PanelMatch(panel.data = dat.anlsys.mmbrs, lag = 2, match.missing = TRUE, 
                            refinement.method = "mahalanobis",
                            covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                              v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                              I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                              e_gdppc + v2paallian, 
                            qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# Propensity score matching
ref.psm.mmbrs <- PanelMatch(panel.data = dat.anlsys.mmbrs, lag = 2, match.missing = TRUE, 
                            refinement.method = "ps.match",
                            covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                              v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                              I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                              e_gdppc + v2paallian, 
                            qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# Propensity score weighing
ref.psw.mmbrs <- PanelMatch(panel.data = dat.anlsys.mmbrs, lag = 2, match.missing = TRUE, 
                            refinement.method = "ps.weight",
                            covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                              v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                              I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                              e_gdppc + v2paallian, 
                            qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)

## compare balance stats
ref.none.mmbrs.blnce <- get_covariate_balance(ref.none.mmbrs, panel.data = dat.anlsys.mmbrs,
                                              covariates = c("v2pavote", "v2elparlel","v2pariglef_ord","v2pawomlab_ord","v2paimmig_ord",
                                                             "e_gdppc", "v2paallian","prmrs.mmbrs.swtch.other","gov.mmbr"))
ref.mal.mmbrs.blnce <- get_covariate_balance(ref.mal.mmbrs, panel.data = dat.anlsys.mmbrs,
                                             covariates = c("v2pavote", "v2elparlel","v2pariglef_ord","v2pawomlab_ord","v2paimmig_ord",
                                                            "e_gdppc", "v2paallian","prmrs.mmbrs.swtch.other","gov.mmbr"))
ref.psm.mmbrs.blnce <- get_covariate_balance(ref.psm.mmbrs, panel.data = dat.anlsys.mmbrs,
                                             covariates = c("v2pavote", "v2elparlel","v2pariglef_ord","v2pawomlab_ord","v2paimmig_ord",
                                                            "e_gdppc", "v2paallian","prmrs.mmbrs.swtch.other","gov.mmbr"))
ref.psw.mmbrs.blnce <- get_covariate_balance(ref.psw.mmbrs, panel.data = dat.anlsys.mmbrs,
                                             covariates = c("v2pavote", "v2elparlel","v2pariglef_ord","v2pawomlab_ord","v2paimmig_ord",
                                                            "e_gdppc", "v2paallian","prmrs.mmbrs.swtch.other","gov.mmbr"))
print(c(ref.none.mmbrs.blnce, ref.mal.mmbrs.blnce, ref.psm.mmbrs.blnce, ref.psw.mmbrs.blnce))

## plot balance before and after refinement (by each method)
# generate dataframe for plot + reshape for plot
ref.none.mmbrs.blnce.df <- as.data.frame(ref.none.mmbrs.blnce)
ref.none.mmbrs.blnce.df <- setDT(ref.none.mmbrs.blnce.df, keep.rownames = TRUE)[]
ref.none.mmbrs.blnce.long <- gather(ref.none.mmbrs.blnce.df, "covariate", "balance",
                                    ref.none.mmbrs.att.v2pavote:ref.none.mmbrs.att.gov.mmbr)
ref.none.mmbrs.blnce.long$rn <- factor(ref.none.mmbrs.blnce.long$rn, levels=c("t_2","t_1","t_0"))
ref.none.mmbrs.blnce.long$dv <- ifelse(
  ref.none.mmbrs.blnce.long$covariate=="ref.none.mmbrs.att.v2pavote",1,0)
ref.none.mmbrs.blnce.long <- filter(ref.none.mmbrs.blnce.long,rn!="t_0")

ref.mal.mmbrs.blnce.df <- as.data.frame(ref.mal.mmbrs.blnce)
ref.mal.mmbrs.blnce.df <- setDT(ref.mal.mmbrs.blnce.df, keep.rownames = TRUE)[]
ref.mal.mmbrs.blnce.long <- gather(ref.mal.mmbrs.blnce.df, "covariate", "balance",
                                   ref.mal.mmbrs.att.v2pavote:ref.mal.mmbrs.att.gov.mmbr)
ref.mal.mmbrs.blnce.long$rn <- factor(ref.mal.mmbrs.blnce.long$rn, levels=c("t_2","t_1","t_0"))
ref.mal.mmbrs.blnce.long$dv <- ifelse(
  ref.mal.mmbrs.blnce.long$covariate=="ref.mal.mmbrs.att.v2pavote",1,0)
ref.mal.mmbrs.blnce.long <- filter(ref.mal.mmbrs.blnce.long,rn!="t_0")


ref.psm.mmbrs.blnce.df <- as.data.frame(ref.psm.mmbrs.blnce)
ref.psm.mmbrs.blnce.df <- setDT(ref.psm.mmbrs.blnce.df, keep.rownames = TRUE)[]
ref.psm.mmbrs.blnce.long <- gather(ref.psm.mmbrs.blnce.df, "covariate", "balance",
                                   ref.psm.mmbrs.att.v2pavote:ref.psm.mmbrs.att.gov.mmbr)
ref.psm.mmbrs.blnce.long$rn <- factor(ref.psm.mmbrs.blnce.long$rn, levels=c("t_2","t_1","t_0"))
ref.psm.mmbrs.blnce.long$dv <- ifelse(
  ref.psm.mmbrs.blnce.long$covariate=="ref.psm.mmbrs.att.v2pavote",1,0)
ref.psm.mmbrs.blnce.long <- filter(ref.psm.mmbrs.blnce.long,rn!="t_0")

ref.psw.mmbrs.blnce.df <- as.data.frame(ref.psw.mmbrs.blnce)
ref.psw.mmbrs.blnce.df <- setDT(ref.psw.mmbrs.blnce.df, keep.rownames = TRUE)[]
ref.psw.mmbrs.blnce.long <- gather(ref.psw.mmbrs.blnce.df, "covariate", "balance",
                                   ref.psw.mmbrs.att.v2pavote:ref.psw.mmbrs.att.gov.mmbr)
ref.psw.mmbrs.blnce.long$rn <- factor(ref.psw.mmbrs.blnce.long$rn, levels=c("t_2","t_1","t_0"))
ref.psw.mmbrs.blnce.long$dv <- ifelse(
  ref.psw.mmbrs.blnce.long$covariate=="ref.psw.mmbrs.att.v2pavote",1,0)
ref.psw.mmbrs.blnce.long <- filter(ref.psw.mmbrs.blnce.long,rn!="t_0")

# plot balance (Figure A1)
ref.none.mmbrs.blnce.plot <- ggplot(dat=ref.none.mmbrs.blnce.long, aes(x=rn, y=balance, alpha=as.factor(dv),
                                                                       group=covariate)) +
  geom_line(aes(size=as.factor(dv))) + 
  geom_hline(yintercept=0.25, linetype="dashed") +
  geom_hline(yintercept=-0.25, linetype="dashed") +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_x_discrete(labels=c("t-2", "t-1")) +
  scale_alpha_discrete(range = c(0.3, 1), labels=c("Covariates","Vote Share")) +
  scale_size_discrete(range = c(0.5,1), labels=c("Covariates","Vote Share")) + 
  ylim(-1.6,1.6) + 
  labs(title="No Refinement", x="Elections Before Treatment", y="Standardzies Mean \nDifferences") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5),
  )

ref.mal.mmbrs.blnce.plot <- ggplot(dat=ref.mal.mmbrs.blnce.long, aes(x=rn, y=balance, alpha=as.factor(dv),
                                                                     group=covariate)) +
  geom_line(aes(size=as.factor(dv))) + 
  geom_hline(yintercept=0.25, linetype="dashed") +
  geom_hline(yintercept=-0.25, linetype="dashed") +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_x_discrete(labels=c("t-2", "t-1")) +
  scale_alpha_discrete(range = c(0.3, 1), labels=c("Covariates","Vote Share")) +
  scale_size_discrete(range = c(0.5,1), labels=c("Covariates","Vote Share")) + 
  ylim(-1.6,1.6) + 
  labs(title="Mahalanobis", x="Elections Before Treatment", y="Standardzies Mean \nDifferences") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5),
  )

ref.psm.mmbrs.blnce.plot <- ggplot(dat=ref.psm.mmbrs.blnce.long, aes(x=rn, y=balance, alpha=as.factor(dv),
                                                                     group=covariate)) +
  geom_line(aes(size=as.factor(dv))) + 
  geom_hline(yintercept=0.25, linetype="dashed") +
  geom_hline(yintercept=-0.25, linetype="dashed") +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_x_discrete(labels=c("t-2", "t-1")) +
  scale_alpha_discrete(range = c(0.3, 1), labels=c("Covariates","Vote Share")) +
  scale_size_discrete(range = c(0.5,1), labels=c("Covariates","Vote Share")) + 
  ylim(-1.6,1.6) + 
  labs(title="Propensity Score\nMatching", x="Elections Before Treatment", y="Standardzies Mean \nDifferences") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5),
  )

ref.psw.mmbrs.blnce.plot <- ggplot(dat=ref.psw.mmbrs.blnce.long, aes(x=rn, y=balance, alpha=as.factor(dv),
                                                                     group=covariate)) +
  geom_line(aes(size=as.factor(dv))) + 
  geom_hline(yintercept=0.25, linetype="dashed") +
  geom_hline(yintercept=-0.25, linetype="dashed") +
  geom_hline(yintercept=0, linetype="dotted") +
  scale_x_discrete(labels=c("t-2", "t-1")) +
  scale_alpha_discrete(range = c(0.3, 1), labels=c("Covariates","Vote Share")) +
  scale_size_discrete(range = c(0.5,1), labels=c("Covariates","Vote Share")) + 
  ylim(-1.6,1.6) + 
  labs(title="Propensity Score\nWeighting", x="Elections Before Treatment", y="Standardzies Mean \nDifferences") +
  theme_bw() +
  theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5),
  )

plot.mmbrs.blnce <- ref.none.mmbrs.blnce.plot + ref.mal.mmbrs.blnce.plot + ref.psm.mmbrs.blnce.plot + ref.psw.mmbrs.blnce.plot +  
  plot_layout(ncol = 2, guides = "collect") & theme(legend.position = 'bottom')
# save plot
ggsave(file="FigureA2.pdf", plot=plot.mmbrs.blnce, width=7, height=7)

#######################################################################
# matched sets - details and size frequency distributions (Figure A3) #
#######################################################################
## calculate matching by selected refinement methods
# moderate inclusiveness measure
PM.results.dlgts <- PanelMatch(panel.data = dat.anlsys.dlgts, lag = 2, match.missing = TRUE, 
                               refinement.method = "ps.weight",
                               covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                 v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                 I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                                 e_gdppc + v2paallian, 
                               qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)
# strict inclusiveness measure
PM.results.mmbrs <- PanelMatch(panel.data = dat.anlsys.mmbrs, lag = 2, match.missing = TRUE, 
                               refinement.method = "ps.weight",
                               covs.formula = ~ I(lag(v2pavote, 1:2)) + 
                                 v2pariglef_ord + v2pawomlab_ord + v2paimmig_ord + v2elparlel +
                                 I(lag(gov.mmbr,1)) + I(lag(prmrs.mmbrs.swtch.other,1)) + 
                                 e_gdppc + v2paallian, 
                               qoi = "att" , lead = 0, forbid.treatment.reversal = FALSE)

## examine treatment units and matched sets
msets.dlgts <- PM.results.dlgts$att    
summary(msets.dlgts)
msets.mmbrs <- PM.results.mmbrs$att    
summary(msets.mmbrs)


# side by side plot (Figure A3)
pdf(file="FigureA3.pdf", width=8, height=3.5)
par(mfrow=c(1,2))
plot(PM.results.dlgts, include.empty.sets = TRUE, 
     main = "Moderate Inclusiveness",
     ylim=c(0,6), xlab = "", ylab = "")
title(xlab = "Matched Set Size", line = 2.5) 
title(ylab = "Frequency", line = 2.5) 
plot(PM.results.mmbrs, include.empty.sets = TRUE,
     main = "Strict Inclusiveness",
     ylim=c(0,6), xlab = "", ylab = "")
title(xlab = "Matched Set Size", line = 2.5) 
title(ylab = "Frequency", line = 2.5) 
dev.off()
